#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use chromosomes;
use linewrap;
use Data::Dump qw(dump);
use db_parameters;
use ortho_species;

print <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    paralogues full functional analysis

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

#
#	Connect to genome database
#
my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);
$dbh->{RaiseError} = 0;
#my @db_param = get_db_parameters();



#_________________________________________________________________________________________

#	get_clusters_above_size


#	Retrieve clusters in order of cluster sizes

#_________________________________________________________________________________________
sub get_closest_other_species($$$\%)
{
	my ($dbh, $species, $protocol_id, $closest_homolog_other_species) = @_;
	my @clusters;
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
            gene_id, 
            other_species_gene_id, 
            cos.ds 
        FROM 
            ortholog_sets JOIN 
            closest_other_species cos USING (gene_id, species, protocol_id) 
        WHERE 
            protocol_id = $protocol_id AND
            ortholog_type ~ '^Orphan inparalog' AND 
            species = '$species';
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for (@$data)
	{
		$closest_homolog_other_species->{$_->[0]} = [$_->[1], $_->[2]];
	}
}




#_________________________________________________________________________________________

#	get_clusters_above_size


#	Retrieve clusters in order of cluster sizes

#_________________________________________________________________________________________
sub get_clusters_above_size($$$)
{
	my ($dbh, $cluster_size, $protocol_id) = @_;
	my @clusters;
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
        ortholog_id  
    FROM 
        orthologs.ortholog_sets 
    WHERE
        ortholog_type ~ 'paralog' AND
        protocol_id = $protocol_id
    GROUP BY ortholog_id
	HAVING count(*) > $cluster_size;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for (@$data)
	{
		push(@clusters, $_->[0]);
	}
	return @clusters;
}


#_________________________________________________________________________________________
#
#	retrieve_interpro_data

#		Retrieve interpro names and descriptions for all matches

#_________________________________________________________________________________________
sub retrieve_interpro_data($\%)
{
	my ($dbh, $interpro_data) = @_;
	print STDERR "\tRetrieve interpro names and descriptions for all matches...\n";
    my $sql_cmd = <<"PL/SQLCMD";
    SELECT
        interpro_ac,
        interpro_name,
        interpro_description
    FROM
        domains.interpro_data;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
    {
        $interpro_data->{$line->[0]} = [$line->[1], $line->[2]];
    }
}




#_________________________________________________________________________________________
#
#	retrieve_other_taxon_prot_interpro_matches

#_________________________________________________________________________________________
sub retrieve_other_taxon_prot_interpro_matches($$\%)
{
	my ($dbh, $protocol_id, $interpro_per_ortholog) = @_;
	print STDERR "\tObtain interpro matches for each ortholog set via homologues in the other species:\n";
	print STDERR "\tSort interpro matches from $ortho_name1, $ortho_name2...\n";
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
            DISTINCT ON (ortholog_id, interpro_ac) 
            ortholog_id, interpro_ac 
        FROM 
            ortholog_sets os JOIN 
            closest_other_species cos USING (protocol_id, gene_id, species) JOIN 
            ens_id e ON (cos.other_species_gene_id = e.gene_id) JOIN 
            taxon.prot_interpro_matches pim ON (e.prot_id = pim.prot_id)
    	WHERE
            ortholog_type ~ '^Orphan inparalog' AND 
            protocol_id =$protocol_id;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		push(@{$interpro_per_ortholog->{$line->[0]}}, $line->[1]);
	}
}


#_________________________________________________________________________________________
#
#	retrieve_taxon_prot_interpro_matches

#_________________________________________________________________________________________
sub retrieve_taxon_prot_interpro_matches($$\%)
{
	my ($dbh, $protocol_id, $interpro_per_ortholog) = @_;
	print STDERR "\tObtain interpro matches for each ortholog set:\n";
	print STDERR "\tSort interpro matches from $ortho_name1, $ortho_name2...\n";
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
		DISTINCT ON (protocol_id, ortholog_id, interpro_ac) 
		ortholog_id, interpro_ac
	FROM
		taxon.prot_interpro_matches NATURAL JOIN
        ens_id JOIN 
		orthologs.ortholog_sets USING (species, gene_id)
	WHERE
        ortholog_type ~ '^Orphan inparalog' AND 
        protocol_id =$protocol_id;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		push(@{$interpro_per_ortholog->{$line->[0]}}, $line->[1]);
	}
}


#_________________________________________________________________________________________
#
#	retrieve_go_matches

#_________________________________________________________________________________________
sub retrieve_go_matches($$\%)
{
	my ($dbh, $protocol_id, $go_per_ortholog) = @_;
	print STDERR "\tObtain GO term matches for each ortholog set:\n";
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
            DISTINCT ON (protocol_id, ortholog_id, GO_id)
            ortholog_id,
            GO_id
        FROM 
            taxon.gene_go_matches gtpg
        NATURAL JOIN 
            (
                SELECT
                    GO_id
                FROM
                    orthologs.significant_go_terms
                WHERE 
                    category = 'duplicate' AND 
                    over_representation = TRUE AND
                    protocol_id = $protocol_id
            ) AS a 
        JOIN
            orthologs.ortholog_sets USING (gene_id)
		WHERE 
    		protocol_id = $protocol_id;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		push(@{$go_per_ortholog->{$line->[0]}}, $line->[1]);
	}
}


#_________________________________________________________________________________________
#
#	retrieve_go_matches_per_geneid

#_________________________________________________________________________________________
sub retrieve_go_matches_per_gene_id($\%)
{
	my ($dbh, $go_per_gene_id) = @_;
	print STDERR "\tObtain GO term matches for each gene_id:\n";
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
            gene_id,
            GO_id
        FROM 
            taxon.gene_go_matches gtpg NATURAL JOIN 
            (
                SELECT
                    GO_id
                FROM
                    orthologs.significant_go_terms
                WHERE 
                    category = 'duplicate' AND 
                    over_representation = TRUE
            ) AS a;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		push(@{$go_per_gene_id->{$line->[0]}}, $line->[1]);
	}
}


#_________________________________________________________________________________________
#
#	retrieve_go_description

#_________________________________________________________________________________________
sub retrieve_go_description($\%)
{
	my ($dbh, $go_descriptions) = @_;
	print STDERR "\tRetrieve GO descriptions for all matches...\n";
    my $sql_cmd = <<"PL/SQLCMD";
    SELECT
        distinct on (go_id)
        go_id,
        description
    FROM
        orthologs.significant_go_terms NATURAL JOIN 
        go_term_descriptions gtd;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
    {
        $go_descriptions->{$line->[0]} = $line->[1];
    }
}






#_________________________________________________________________________________________

#	print_go_data

#_________________________________________________________________________________________
sub print_go_data($$$\%\%\%$)
{
	my ($dbh, $ortholog_id, $printH, $go_descriptions, $go_per_ortholog, $overrepresented_GOids, $species) = @_;
	return unless $go_per_ortholog->{$ortholog_id};
	my @matches = @{$go_per_ortholog->{$ortholog_id}};
	return unless @matches;

	print $printH "    GO matches for this cluster:\n";
	for my $go_id(@matches)
	{
		my $over = "";
		if (exists $overrepresented_GOids->{$go_id}{$species})
		{
			$over = ($overrepresented_GOids->{$go_id}{$species}[1]) ? "++over++" : "__over__";
		}
		die "[$go_id]\n" unless exists $go_descriptions->{$go_id};
		die "[$go_id]\n" unless defined $go_id;
		die "[$go_id]\n" unless defined $go_descriptions->{$go_id};
		printf  $printH "      %-15s %8s %-20s\n", "GO:$go_id", $over, $go_descriptions->{$go_id};
	}
	print  $printH  "\n";
}


#_________________________________________________________________________________________
#
#
#	retrieve_overrepresented_GO_ids
#
#_________________________________________________________________________________________
sub retrieve_overrepresented_GO_ids($$\%)
{
	print STDERR "\tGet significantly overrepresented GO IDs...\n";
	my ($dbh, $protocol_id, $overrepresented_GO) = @_;
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		GO_id,
        species,
        CASE WHEN go_slim THEN 1 ELSE 0 END 
	FROM 
		orthologs.significant_go_terms
    WHERE
        category = 'duplicate' AND 
        significance < 0.0001 AND 
        over_representation = TRUE AND
        protocol_id = $protocol_id
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";

	# save GO identifiers per GO category per gene id
	for my $line (@$data)
	{
		++$overrepresented_GO->{$line->[0]}{$line->[1]}[$line->[2]];
	}

}



		

#_________________________________________________________________________________________

#	print_interpro_data

#_________________________________________________________________________________________
sub print_interpro_data($$$\%\%)
{
	my ($dbh, $ortholog_id, $printH, $interpro_data, $interpro_per_ortholog) = @_;
	return unless $interpro_per_ortholog->{$ortholog_id};
	my @matches = @{$interpro_per_ortholog->{$ortholog_id}};
	return 0 unless @matches;

	print $printH "    Interpro matches for this cluster:\n";
	for my $interpro_ac(@matches)
	{
		printf  $printH "      %-15s %-20s %s\n", $interpro_ac, @{$interpro_data->{$interpro_ac}};
	}
	print  $printH  "\n";
	return 1;
}





#_________________________________________________________________________________________

#	gene loci per cluster

#		clusters are grouped by ortholog_id, species and adjacent clusters

#_________________________________________________________________________________________
sub get_per_cluster_loci($\%$$)
{
	my ($dbh, $cluster_loci, $species, $protocol_id) = @_;
	print STDERR ("\tRetrieving clusters for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		  ortholog_id,
		  local_cluster_id,
		  gene_id,
		  prot_id,
		  gene_position,
		  chromosome
	FROM
        orthologs.ortholog_sets NATURAL JOIN
		orthologs.local_paralogs NATURAL JOIN
		taxon.gene_loci NATURAL JOIN 
        orthologs.real_genes
    WHERE
        protocol_id = $protocol_id AND 
        ortholog_type ~ 'paralog' AND
        species = '$species'
	ORDER BY ortholog_id, local_cluster_id, gene_position;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		# ortholog_id, species, local_cluster_id, [gene_id, prot_id, gene_position, chromosome]
		push(@{$cluster_loci->{$line->[0]}{$species}{$line->[1]}}, [$line->[2],$line->[3],$line->[4],$line->[5]]);
	}
}

#_________________________________________________________________________________________

#	get_ortholog_cluster_ks

#_________________________________________________________________________________________
sub get_ortholog_cluster_ks($\%$)
{
	my ($dbh, $cluster_ks, $protocol_id) = @_;
	print STDERR ("\tRetrieving dS for each cluster...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
        DISTINCT ON (ortholog_id)
        ortholog_id,
        dS
    FROM
    	orthologs.ortholog_sets
    WHERE
        ortholog_type ~ 'paralog' AND
        protocol_id = $protocol_id 
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$cluster_ks->{$line->[0]} =$line->[1];
	}
}



#_________________________________________________________________________________________

#	gene description

#_________________________________________________________________________________________
sub get_gene_descriptions($\%$)
{
	my ($dbh, $gene_descriptions, $species) = @_;
	print STDERR ("\tRetrieving gene descriptions for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		  gene_id,
		  description
    	FROM
    		taxon.gene_descriptions NATURAL join
            ens_id
        WHERE 
            species = '$species';    
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$gene_descriptions->{$line->[0]} =$line->[1];
	}
}

#_________________________________________________________________________________________

#	get_ortholog_gene_descriptions
#		treat genes also like ortholog sets

#_________________________________________________________________________________________
sub get_ortholog_gene_descriptions($\%$)
{
	my ($dbh, $gene_descriptions, $species) = @_;
	print STDERR ("\tRetrieving gene descriptions for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
            DISTINCT ON (ortholog_id, description)
            ortholog_id,
            description
    	FROM
    		taxon.gene_descriptions NATURAL join
            ortholog_sets
        WHERE 
            species = '$species';
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		push(@{$gene_descriptions->{$line->[0]}}, $line->[1]);
	}

	$sql_cmd = <<"PL/SQLCMD";
    SELECT
		  gene_id,
		  description
    	FROM
    		taxon.gene_descriptions NATURAL join
            ens_id
        WHERE 
            species = '$species';    
PL/SQLCMD
	$data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		push(@{$gene_descriptions->{$line->[0]}}, $line->[1]);
	}
}



#_________________________________________________________________________________________

#	protein exon count

#_________________________________________________________________________________________
sub get_prot_exon_counts($\%$)
{
	my ($dbh, $exon_counts, $species) = @_;
	print STDERR ("\tRetrieving prot exon_counts for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
    		  prot_id,
    		  real_exon_count,
              disruptions_count
    	FROM
    		taxon.prot_exon_stats NATURAL join
            ens_id
        WHERE 
            species = '$species';
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$exon_counts->{$line->[0]} =[$line->[1], ($line->[2] ? $line->[2] : '')];
	}
}


#_________________________________________________________________________________________

#	protein exon count

#_________________________________________________________________________________________
sub get_prot_seqlen($\%)
{
	my ($dbh, $prot_len) = @_;
	print STDERR ("\tRetrieving prot lengths (residues) ...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		  prot_id,
		  len
	FROM
		orthologs.ortholog_sets NATURAL JOIN 
        accessions NATURAL join
        sequences;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$prot_len->{$line->[0]} = $line->[1];
	}
}


#_________________________________________________________________________________________

#	gene exon counts

#_________________________________________________________________________________________
sub get_gene_exon_counts($\%$)
{
	my ($dbh, $exon_counts, $species) = @_;
	print STDERR ("\tRetrieving gene exon_counts for $species...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
            gene_id,
            exon_count
        FROM
            taxon.gene_coding_exon_counts NATURAL JOIN 
            ens_id
        WHERE 
            species = '$species';
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$exon_counts->{$line->[0]} =$line->[1];
	}
}




#_________________________________________________________________________________________
#
#	get_cluster_sizes

#	Retrieve clusters in order of cluster sizes

#_________________________________________________________________________________________
sub get_cluster_sizes($\%\@$$)
{
	my ($dbh, $cluster_sizes, $clusters, $species, $protocol_id) = @_;
	print STDERR ("\tRetrieving the list of clusters in descending order of $species cluster size...\n");

	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
        ortholog_id, 
        sp.paralog_count
    FROM
        orthologs.paralog_stats sp NATURAL JOIN 
        (
        SELECT 
                ortholog_id 
            FROM 
                ortholog_sets
            WHERE 
                ortholog_type ~ '^Orphan inparalog'
            GROUP BY 
                ortholog_id
        ) AS a
    WHERE
		sp.paralog_count >= 1 AND
        sp.species = '$species' AND
        protocol_id = $protocol_id
    ORDER BY 
        sp.paralog_count desc;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd. "!";
	for my $line (@$data)
	{
		$cluster_sizes->{$line->[0]} = $line->[1];
		push(@$clusters, $line->[0]);
	}
}

#_________________________________________________________________________________________
#
#	print cluster sizes

#_________________________________________________________________________________________
sub print_cluster_sizes(*\@\%\%$)
{
	my ($printH, $ortholog_ids, $cluster_sizes, $cluster_ks, $species) = @_;
	print $printH "Cluster sizes for clusters >= 4:\n";
	printf $printH "      %10s %5s\n", " size", "dS";
	printf $printH "      %10s %5s\n", "-----", "-----";
	for my $ortholog_id(@$ortholog_ids)
	{
		printf $printH "      %10s %5s\n",
					$cluster_sizes->{$ortholog_id},
					0#$cluster_ks->{$ortholog_id}
			if ($cluster_sizes->{$ortholog_id} >= 4);
	}
	print $printH "\n\n\n";

}



#_________________________________________________________________________________________

#	print all summary info

#_________________________________________________________________________________________
sub print_summary(**\@\%\%\%\%\%\%\%\%\%\%\%\%$\%\%)
{
	my ($printIDS,  $printH,  $ortholog_ids, $cluster_sizes, $exon_counts, $gene_descriptions, 
		$cluster_loci, $cluster_ks, 
		$go_descriptions, $go_per_ortholog,
		$overrepresented_GOids,
		$interpro_data, $interpro_per_ortholog, $go_per_gene_id, $prot_len, $species,
		$closest_homolog_other_species, $other_interpro_per_ortholog) = @_;
	
	my $file_num = 0;
	for my $ortholog_id(@$ortholog_ids)
	{
		last if ($cluster_sizes->{$ortholog_id} < 3);
		my $ortholog_id_printed = 0;
		++$file_num;
		print $printH "Cluster: #$file_num        (ortholog_id = $ortholog_id in panda)\n";
		printf $printH "    %s $species genes dS = %s\n",
								$cluster_sizes->{$ortholog_id},
								0;#$cluster_ks->{$ortholog_id};
		if (!print_interpro_data($dbh, $ortholog_id, $printH, %$interpro_data, 
							%$interpro_per_ortholog))
		{
			print $printH "    Interpro matches via homologues in the other species.\n";
			print_interpro_data($dbh, $ortholog_id, $printH, %$interpro_data, 
								%$other_interpro_per_ortholog);
				#print STDERR dump %$other_interpro_per_ortholog;

		}
		print_go_data($dbh, $ortholog_id, $printH, %$go_descriptions, %$go_per_ortholog, %$overrepresented_GOids, $species);
		printf $printH "      %-22s %8s %5s %5s %4s %6s\n", "gene identifier", "order", "chrm", "exons", "stop", "length";
		printf $printH "      %-22s %8s %5s %5s %4s %6s\n", "-" x 20, "-" x 5, "-" x 4, "-" x 5, "-" x 4, "-" x 6;

		{
			my $prev_gene_id = '';
			for my $sub_cluster_id(sort{$a <=> $b} keys %{$cluster_loci->{$ortholog_id}{$species}})
			{
				for (@{$cluster_loci->{$ortholog_id}{$species}{$sub_cluster_id}})
				{
					my ($gene_id, $prot_id, $gene_pos, $chromo) =  @$_;
					my ($prot_exon_count, $gene_exon_count) = ($exon_counts->{$prot_id}, $exon_counts->{$gene_id});
					if ($gene_id ne $prev_gene_id)
					{
						my @gene_id = split /&/, $gene_id;
						if (@gene_id > 1)
						{
							printf $printH "      %-22s %8d %5s %4s               ", 
										$gene_id[0].'_split', $gene_pos, $chromo, "";
							if (!$ortholog_id_printed)
							{
								print $printIDS "$ortholog_id\t$gene_id[0]\t??\n";
							}
							$ortholog_id_printed = 1;
						}
						else
						{
							printf $printH "      %-22s %8d %5s %4s               ", 
										$gene_id, $gene_pos, $chromo, "";#$gene_exon_count;
							if (!$ortholog_id_printed)
							{
								print $printIDS "$ortholog_id\t$gene_id\t??\n";
							}
							$ortholog_id_printed = 1;
						}

						my $closest_homologue;
						if (exists $closest_homolog_other_species->{$gene_id})
						{
							my $ds;
							($closest_homologue,  $ds) = @{$closest_homolog_other_species->{$gene_id}};
							print $printH "[ds = $ds / $closest_homologue] ";
						}
						if (exists $gene_descriptions->{$gene_id})
						{
							my $desc = $gene_descriptions->{$gene_id};
							$desc =~ s/\[Source:.+\]//;
							print $printH $desc;
						}
						else
						{
							if (defined $closest_homologue &&
								exists $gene_descriptions->{$closest_homologue})
							{
								my $desc = $gene_descriptions->{$closest_homologue};
								$desc =~ s/\[Source:.+\]//;
								print $printH "===> $desc";
							}
						}
						$prev_gene_id =$gene_id;
						print $printH "\n";
					}
					if (exists $exon_counts->{$prot_id})
					{
						my $go_ids = '';
						#$go_ids = join(", ",  map {"GO:".$_} @{$go_per_gene_id->{$gene_id}}) 
						#					if (exists $go_per_gene_id->{$gene_id});
						printf $printH "          %-25s         %5s %4s %6s            %s\n", 
							$prot_id, $prot_exon_count->[0], $prot_exon_count->[1], $prot_len->{$prot_id}, $go_ids;
					}
				}
				print $printH "\n";
			}
		}
		print $printH "\n\n";
	}
	print STDERR "\tFinished summary information...\n";
}













#_________________________________________________________________________________________

#	get data

#_________________________________________________________________________________________

system ("mkdir -p $dir_orthologs_output/paralog_function");


# exon counts
my %exon_counts;
get_prot_exon_counts($dbh, %exon_counts, $ortho_name2);
get_prot_exon_counts($dbh, %exon_counts, $ortho_name1);
get_gene_exon_counts($dbh, %exon_counts, $ortho_name2);
get_gene_exon_counts($dbh, %exon_counts, $ortho_name1);

# gene descriptions
my %gene_descriptions;
get_gene_descriptions($dbh, %gene_descriptions, $ortho_name2);
get_gene_descriptions($dbh, %gene_descriptions, $ortho_name1);






# go data
my %go_descriptions;
retrieve_go_description($dbh, %go_descriptions);
my %go_per_gene_id;
retrieve_go_matches_per_gene_id($dbh, %go_per_gene_id);


# protein length data
my %prot_len;
get_prot_seqlen($dbh, %prot_len);

	
for my $protocol (get_orthologs_protocols())
{
	my ($protocol_name,  $protocol_id) = @$protocol;
	next unless $protocol_id = $curr_protocol_id;
	print STDERR "    Orthology protocol [$protocol_name]...\n";

	#	Overrepresented GO identifiers for each protocol
	#
	my %overrepresented_GOids;
	retrieve_overrepresented_GO_ids($dbh, $protocol_id, %overrepresented_GOids);
	
	
	# go data
	my %go_per_ortholog;
	retrieve_go_matches($dbh, $protocol_id, %go_per_ortholog);
	
	# interpro data
	my (%interpro_data, %interpro_per_ortholog);
	retrieve_interpro_data($dbh, %interpro_data);
	retrieve_taxon_prot_interpro_matches($dbh, $protocol_id, 
										 %interpro_per_ortholog);

	my %other_interpro_per_ortholog;
	retrieve_other_taxon_prot_interpro_matches(	$dbh, $protocol_id, 
												%other_interpro_per_ortholog);

				
	for my $ortho_name($ortho_name1, $ortho_name2)
	{
		my %cluster_loci;
		get_per_cluster_loci($dbh, %cluster_loci, $ortho_name, $protocol_id);

		my %cluster_ks;
		get_ortholog_cluster_ks($dbh, %cluster_ks, $protocol_id);

		my %closest_homolog_other_species;
		get_closest_other_species($dbh, $ortho_name, $protocol_id, %closest_homolog_other_species);



		#
		#	get cluster sizes
		#
		my %cluster_sizes;
		my @cluster_ids;
		get_cluster_sizes($dbh, %cluster_sizes, @cluster_ids, $ortho_name, $protocol_id);

		open ORTHO_IDS,			">$dir_orthologs_output/paralog_function/orphans.$ortho_name.ortholog_ids" or die;
		open SPECIES_RESULTS,	">$dir_orthologs_output/paralog_function/orphans.$ortho_name.function" or die;
		print_cluster_sizes(*SPECIES_RESULTS, @cluster_ids, %cluster_sizes, %cluster_ks, $ortho_name);

		print STDERR "\tPrinting summary information ...\n";
		print_summary(*ORTHO_IDS, *SPECIES_RESULTS, @cluster_ids, %cluster_sizes, %exon_counts, 
						%gene_descriptions, %cluster_loci, %cluster_ks, 
						%go_descriptions, %go_per_ortholog,
						%overrepresented_GOids,
						%interpro_data, %interpro_per_ortholog, 
						%go_per_gene_id, %prot_len, $ortho_name, %closest_homolog_other_species,
						%other_interpro_per_ortholog);
	}
	
}


log_pipeline_stage($dbh, $curr_protocol_id, 46, 'table analyse paralogues', $start_time);
$dbh->disconnect();

